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Abstract 

Recently, the intrinsic sampling method has been developed in order to obtain, from molecular 
simulations, the intrinsic structure of the liquid-vapor interface that is presupposed in the classical 
capillary wave theory. Our purpose here is to study dynamical processes at the liquid-vapor 
interface, since this method allows tracking down and analyzing the movement of surface molecules, 
thus providing, with great accuracy, dynamical information on molecules that are "at" the interface. 
We present results for the coefficients for diffusion parallel and perpendicular to the liquid-vapor 
interface of the Lennard- Jones fluid, as well as other time and length parameters that characterize 
the diffusion process in this system. We also obtain statistics of permanence and residence time. 
The generality of our results is tested by varying the system size and the temperature; for the later 
case, an existing model for alkali metals is also considered. Our main conclusion is that, even if 
diffusion coefficients can still be computed, the turnover processes, by which molecules enter and 
leave the intrinsic surface, are as important as diffusion. For example, the typical time required 
for a molecule to traverse a molecular diameter is very similar to its residence time at the surface. 
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I. INTRODUCTION 



Inhomogeneous systems present a number of features that make them intrinsically more 
complicated than bulk systems. The fact that the equilibrium state of the system depends 
on the position causes a number of physical quantities to be likewise dependent on the 
position (such as the molecular number density), or even ill-defined (such as the pressure 
tensor). This also applies to dynamical properties — the most important one of these, (self- 
diffusion, is complicated by the fact that tracer molecules cannot be followed at will for any 
given length of time, since they will enter and abandon zones that have different dynamical 
properties. 

It is therefore not possible in general to obtain a value for the diffusion coefficient D by the 
well known Einstein relation for the mean standard deviation (MSD) of the displacements 

(r 2 ) > 6Dt. (1) 

t^oo 

The direct application of such a formalism to an inhomogeneous fluid would result in an 
average result containing contributions from the bulk phases and the interfaces. Indeed, D 
can be expected to have different values at different parts of the system. (The same problem 
would arise of course in the other main approach to D: by means of the Green-Kubo formula 
involving the velocity autocorrelation function.) 

This is true in particular for the best-known inhomogeneous fluid system: the liquid- vapor 
interface. In this case, a liquid phase and its vapor are separated by an interfacial region 
which, on average, is flat. The spacial dependence is therefore limited in this case to one 
Cartesian coordinate, which we will take as z. For example, the mean density profile, p(z), 
is obtained from simulation by defining slab in the z direction ("binning"), and collecting 
occupation statistics for each of the slabs. This way, a profile is obtained that shows two 
plateaus at constant values corresponding to the liquid and vapor densities, and a typically 
monotonic interfacial variation between them. Theoretical approaches, from the pioneering 
van der Waals theory to the most recent density functional approximations, may be used 
to directly obtain p(z), which depends only on the temperature, T, and on the molecular 
interactions. The density profiles may be much more structured in other, apparently more 
complex, systems like a dense fluid against a planar wall potential, but the apparent sim- 
plicity of the liquid-vapor interface hides a much deeper difficulty.-^ The fluctuations of a 
free liquid surface have capillary wave (CW) modes with very low frequencies, and hence 
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low excitation energies for long wavelengths. In the absence of any external potential, the 
thermodynamic limit of a macroscopic free liquid surface becomes undetermined, and the 
interface would be fully delocalized by the long wavelength CWs. The Earth gravity field, 
which would be fully irrelevant for the thermodynamic properties of one-phase systems up 
to the scale of meters, becomes crucial to stabilize the liquid surface, damping the CW 
fluctuations for wavelengths larger than millimeters, and amplitudes larger than about one 
molecular diameter. Still, the mean density profile under Earth gravity conditions would be 
smoother than the one observed with typical computer simulations, which employ transverse 
box sizes in the range of 10 — 30 molecular diameters (the interfacial width can be estimated 
to be about twice as large in the first case 3 ). Within that limited range for the transverse size 
the inclusion of the Earth gravity would be irrelevant, but it is already possible to observe 
the dependence of p(z) with the transverse size of the simulation box.-^ 

The CW fluctuations are a severe nuisance for the study of the molecular diffusion at 
the liquid surface, since the molecular kinetics, relevant for any physicochemical process 
at the surface, is mixed with large collective fluctuations for the instantaneous position of 
the surface, which give a smooth and size dependent p{z). The simplest estimations of 
surface diffusion properties consider the molecules within thin slabs placed in the region 
with inhomogeneous values of p(z),—£ in some cases separating diffusion in parallel and 
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on liquid- vapor interfaces, and Refs. 
liquids (typically, water) adsorbed on, or confined in, different substrates. Another approach 
is to consider some operational definition of the outmost liquid molecules and track their 
dynamics for a limited time.— Such procedures can only yield a coarse distinction between 
bulk and surface properties, given the obvious arbitrariness in the choice of parameters such 
as the surface slab width, the definition of outmost molecules, the complications associated 
with the diffusing molecules leaving and entering the different domains, and the blurring 
effect of the area-dependent CW fluctuations. 

Classical capillary wave theory (CWT) 1 gives a framework to interprete the effects of 
the CW fluctuations, and provides an accurate extrapolation of the p(z) profile obtained in 
typical computer simulations to larger sampling sizes, including the effects of weak gravity 
fields. However, only over the last decade has CWT become a practical tool to extract the 
intrinsic molecular properties of a liquid surface, from the broad distributions produced by 
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the CW fluctuations. The theory assumes that an intrinsic surface (IS) may be defined, 
to describe the instantaneous boundary between the two coexisting phases, so that the 
molecular distribution referred to that surface would give an intrinsic density profile sharper 
that p(z) and, more importantly, independent of the transverse sampling size. Over the 
last decade, the increasing resolution in X-ray reflectivity data allowed the deconvolution 
of the Gaussian CW distribution out of the surface structure factor, and hence to obtain 
experimental results for the intrinsic profile in cold liquid metal surfaces, with a clear atomic 
layering structure.— Similar results were obtained in computer simulations of simple fluid 
models,— 1 ^ whenever the frustration of the freezing allowed to explore low temperatures, 
T/T c < 0.2, relative to the critical one. As in experiments, the deconvolution is possible 
because even p(z) presents some layering at these low temperatures. These results indicate 
that the typical smooth shape of p(z) in a LV interface results from the convolution of the 
Gaussian CW fluctuations with a strongly layered intrinsic profile, such that it should be 
possible to identify, with reasonable confidence, the outmost molecular layer of the liquid 
phase. More recently, that concept led to the development of intrinsic sampling method in 
computer simulations,—^ 1 ^ based on the operational definition of the CWT intrinsic surface 
as a geometrical locus for that first molecular layer, so that the intrinsic profile, and other 
molecular intrinsic properties of the liquid surface,— 1 ^ may be extracted with accuracy and 
high resolution, without the CW blurring observed in p(z). The intrinsic profile calculated 
with such methods represent a direct, quantitative link between the generic framework of 
the CWT and the computer simulations of liquid surfaces. The result is a deeper structural 
understanding of the surface: the capillary wave fluctuations, which cause an area-dependent 
broadening of the mean profile p(z), are absent in the intrinsic one. These two profiles are, 
in general, remarkably different: between the plateaus corresponding to the liquid and vapor 
densities the intrinsic one is far from monotonic, showing a marked layering. In fact, it much 
more resembles the pair correlation function or the density profiles close to a hard wall. 

The aim of our work here is to explore the application of the intrinsic sampling method to 
the analysis of the molecular diffusion at liquid surfaces. We may track down and analyze the 
movements of the surface molecules, i.e. those that define the IS. Thus we may obtain, with 
great accuracy, dynamical information on molecules that are "at" the interface, without the 
arbitrarity in the choice of a surface slab, and making it possible to separate the molecular 
diffusion on the surface from the fluctuations on the local position of the surface. The reader 
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is referred to the previous references for a complete description of the method, its variations, 
and their results. The only relevant aspect here is that for each instantaneous configuration 
of the system, this method selects a set of surface molecules, i.e. those identified as belonging 
to the outmost liquid layer, and called the IS pivots in the above references. Within the 
intrinsic disorder of a liquid surface we cannot expect that such molecular layer would have 
a sharp definition — it should indeed be regarded as a soft but robust concept: changes in 
the set of parameters used for the operational IS definition would produce small changes 
in the set of molecules which are identifyed as belonging to the surface, but (at least for 
T/T c < 0.8) there is a clear optimal choice within rather narrow windows for all parameters 
that must be fine tuned, as commented at the end of the next section. In particular, the 
two-dimensional density of the first liquid layer, i.e. the number of surface molecules per 
unit area, is fairly well defined and provides useful information on the molecular structure 
of a liquid surface, which is blurred in the usual description in terms of p(z)M 

We begin with a section on methodology, Section [Til on both the simulation details and 
the way the IS is calculated. We then consider diffusion in Section IHIl divided in three 
subsections that go from the best known case, bulk diffusion, through diffusion parallel 
to the interface, to the more involved case of diffusion perpendicular to it. We further 
analyze our results in Section IIVI where we consider systems at other temperatures and 
other transverse areas. We conclude with some remarks in [V] 

II. METHODOLOGY 



We consider the standard Lennard- Jones fluid, in which molecules interact through a 
pairwise potential of the form 




In order to obtain dynamical information, molecular dynamics (MD) simulations are 
performed, using the software package dl_poly.— The systems consists of a thick slab of 
liquid surrounded by vapor. We therefore obtain two LV interfaces. Since the two are 
independent (if the liquid is thick enough), the properties measured in both should be 
averaged in order to improve the accuracy, but for the sake of clarity we will just present 
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results obtained in one of the LV interfaces. 

Our systems consist of 2592 molecules, generally (unless otherwise indicated) at a tem- 
perature of k^T = 0.678e, which is our estimate for the triple point temperature of the LJ 
fluid (at this cutoff, see Ref. |42| for a discussion of the effect of truncation on this tempera- 
ture). The simulation cell is a square box of dimensions L x L x L z , with (unless otherwise 
indicated) L = 10.46a and L z = 90a. Periodic boundary conditions are employed in all 
three directions. The time step is set at a reduced value of dt = 4.56 x 10~ 3 o~A/m/e. The 
systems starts from either a crystalline configuration or a system at other temperature and 
are equilibrated for 10 6 time steps in the NVT ensemble (with a Nose-Hoover thermostat 
with a time constant lOdt). 

After this period, configurations are obtained in the NVE ensemble, in order to eliminate 
any possible spurious effects of the thermostat on the dynamics (in any case, we have checked 
these effects are typically negligible). For the results presented here, 50000 configurations 
are analyzed. We have found that only one configuration out of ten needs to be analyzed, 
since the dynamics is still slow for a time step of 10<it. The analysis of the MSD is carried out 
in the standard way; more sophisticated treatments aimed at reducing computer storag ei2 
are not needed in this case. 



For each of these configurations, an IS analysis is carried out, as described in Ref. [37, 
requiring two operational parameters. One of them, n s , has a very clear physical significance: 
the two-dimensional density of the IS, i.e. the number of surface molecules to be selected 
per unit area. Here we set n s pa 0.8a~ 2 , with a twofold motivation; on one hand, this value 
had already been identified as the most physical one from close inspection of the intrinsic 
profile for the LJ fluid.— On the other, the same set of simulations described here can be 
used to independently obtain a value for this parameter, confirming this value as the best 
one from a kinetic analysis (see Ref.Q for details on this procedure). The other parameter, 
q n , fixes the maximum wave number to be used in the Fourier representation of the IS. 
This mathematical surface associated to each molecular configuration within the intrinsic 
sampling method will not appear explicitly this article, but still it is essential for the self- 
consistent procedure used to select the IS molecules. We have combined here a basis of 12 
planar waves for each x and y direction, which with the transverse box size L = 10.46a 
means a value of q n = 1.15 x 27r/a, within the optimal range of choices as explained in the 
above references. 
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FIG. 1: Mean fraction of molecules which remain at the intrinsic surface for, at least, a time 
t, N s (t)/N s (0), versus reduced time t, (units of cr-v/m/e). Dotted line: linear regression of the 
exponential decay at later times, shown to intercept t = at about 0.8. 

In each analyzed configuration we select a number of surface molecules, N s = n s L 2 
(typically 88 is our optimal choice here) which self-consistently define the IS. The averages 
of the surface self-diffusion properties require us to follow the individual history of each of 
these molecules, i.e. to identify their permanence, exit, and possible reentrance in the IS list. 
All the reported surface diffusion properties are associated to the movement of molecules 
that stay continuously in the list of IS molecules, and this requirement sets a clear limitation 
on the available sampling times. In Figure HJ N s (t) represents the mean number of molecules 
which remain at the IS for, at least, a time t; the exponential decay N s (t) ~ exp(— t/r), 
with a typical decay time, or residence time of r m 4.3a \frnft (at k-^T = 0.678e), which we 
also list in Table [H This treatment mirrors the classical definition of the bulk residence time 
as the exponential decay time for the process by which neighboring molecules drift apart.— 
It is also useful to keep in mind that for the case of Argon with the usual LJ parameters 
c/&b = 119. 8K and a ~ 0.3405nm one unit of reduced time corresponds to 2.15 picoseconds. 
Thus, the decay time would be approximately 9.2ps for Argon. 

This short value for the residence time implies that the sampling size would rapidly 
decrease (and the statistical noise increase) for large t/r. The extrapolation to t = of this 
exponential decay at longer times indicates that about 80% of the molecules at the IS will 
leave it at the constant rate of 1/r, associated to the long time exponential decay, while 
the remaining 20% of the molecules get out of the IS much faster, sometimes to undergo 
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FIG. 2: Mean density profile, p(z) (solid line), density profile of surface molecules, p s {z) (dotted 
line), and density profile of old surface molecules, uqP {z) (dashed line). The profiles are given in 
reduced units of <7~ 3 . 

a rapid reentrance. These molecules may be regarded as those which are only marginally 
associated to the outmost liquid layer, e.g. those which may be interpreted either as a 
local intrusion of the IS towards the bulk liquid, or alternatively as a local extrusion of the 
next-outmost layer towards the surface. The disordered structure of the liquid makes the 
existence of such ambiguities unavoidable, and different recipes to identify the IS from the 
molecular positions could make different assignments to those molecules to be in, or out 
of, the list of surface molecules. That is what makes the concept of the outmost liquid 
layer a soft one, but at the same time a rather robust one, since any reasonable choice of 
the tunable parameters would agree in the selection of the large majority of the surface 
molecules. For the diffusion properties analyzed here there is a further advantage, since 
the relevant information to get the effective surface diffusion coefficients comes from the 
sampling of molecular displacements at the longest possible times, which automatically 
selects the properties of those molecules with long permanence at the surface, representing 
the outmost liquid layer with any sensible choice for the parameters in the IS definition. 
The only practical inconvenience of the rapid turnover of some surface molecules is that the 
time interval between analyzed configurations, At, has to be short enough to make sure that 
the estimation of N s (t) is not affected by further reduction — in practice that is achieved, 
at k^T = 0.678e, with At = lOdt, when only about two and a half surface molecules are 
changed on average between consecutive configurations.— 
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The mean profile of the whole system, p(z), and that restricted to the surface molecules, 
p s (z), are presented in Figure [2], the later having a typical Gaussian shape which represents 
the fluctuations of the IS, and which becomes wider with increasing temperature and trans- 
verse sizes of the simulation box. Notice that despite the mean profile character of p s (z), 
all the selected surface molecules lie exactly on the instantaneous IS, we may therefore use 
them to follow exactly the self- diffusion of the molecules at the outmost liquid layer, without 
the coarsening effect of the CW fluctuations if we were selecting the molecules within a fixed 
surface slab. Moreover, we may keep track of the mean profiles for surface molecules which 
have continuous permanence in the IS list for more than a time t, thus defining a distribu- 
tion p s (z, t), whose integral over z is directly linked to the number of surface molecules older 
than t, 



We may expect that any diffusion property sampled for relatively large times would corre- 
spond to molecules distributed as 



where the time dependences factorizes into the same exponential decay as N s (t), and we 
define a ^-distribution of old surface molecules, P {z), which is normalized to unity. The 
prefactor has been discussed above to be n « 0.8n s , i.e. representing about 80% of the 
surface molecules. This is indeed the case, and in Figure [2] we compare the mean profile 
of the whole set of surface molecules, p B (z) = p 8 (z, 0), and the distribution of old surface 
molecules normalized to n , i.e., n P Q (z). The distribution of old surface molecules is 
narrower, and asymmetric with respect to the whole one, and that may be interpreted in 
terms of the rate at with molecules are incorporated to, or deleted from, the list of surface 
molecules, as a function of their position z. 

Since we select a fixed number of surface molecules in each configuration, the loss of 
molecules from the IS layer, either towards the liquid or the vapor sides, is always com- 
pensated by the incorporation of new ones, and the time reversal symmetry of the MD 
guarantees that the inflow and outflow of surface molecules at a given value of z are iden- 
tical. We denote by U{ (z) that input/output rate per molecule, which may be directly 
sampled along our simulations and are presented in Figure [3], together with P (z), for ease 
of comparison. The results are fairly well compatible with rates being independent of the 




(3) 



p s (z,t) = n P o (z)exp(-t/T) 



(4) 
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FIG. 3: Distribution of the input/output rate per molecule, U[ (z), normalized to the proper value 



set by Equation [6] (hence, reduced in units of inverse time, \J e/m/cr). (solid line), together with 
P (z) in reduced units of a~ 1 , normalized to unity (dashed line). 

previous permanence time of the molecule at the IS, so that 

d f 

-iV s (t) = —L 2 / dz Ps {z,t)u io {z), (5) 



which for large t implies 



1 1 dN s (t) 



dzP Q (z)u io (z). (6) 



r N s (t) dt 

The shape of v- lo {z) is clearly composed of two contributions: a large peak corresponding 
to the exit /entrance of surface molecules to/from the bulk liquid, an a much smaller one 
corresponding to the exit /entrance of surface molecules to/from the bulk vapor. In the 
following section we show how to extract the surface diffusion coefficients from the informa- 
tion from the observed displacements of the surface molecules, and the information on their 
turnover distributions contained in P Q (z) and Vi Q {z). 



III. DIFFUSION IN THE BULK, AND AT THE INTRINSIC SURFACE 

We discuss diffusion for the particular choice of temperature ksT = 0.678e, in three 
different cases: bulk diffusion, diffusion parallel to the interface, and diffusion perpendicular 
to it. 
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FIG. 4: Square root of the mean standard deviation of displacements for different subset of 



molecules versus time, in reduced units: a for the square root of MSDs and o^J m/e for time. 
Solid lines: x, y, and z components in the bulk liquid (hardly distinguishable), long-dashed lines: 
x and y components for the intrinsic surface (hardly distinguishable), short-dashed line: z compo- 
nent for the intrinsic surface. Dotted lines: prediction from the Maxwellian velocity distribution, 
Eq. (0), and fits to Einsteinian diffusion equations, Eq. Q, with two values of the self diffusion 
coefficient, one corresponding to the bulk liquid, the other to the intrinsic surface. 

A. Diffusion in the bulk 

In Figure H] we plot the square root of the MSD for the x, y, and z components of the 
displacement versus time for the bulk liquid phase. At short times all curves tend to the 
same line, with a slope of 2. This is the ballistic regime: times so short that collisions can 
be neglected. In this case the Maxwellian distribution for velocities directly provides a value 
for the MSD 

(x 2 ) = ^t 2 , (7) 

which is the dotted line in the graph. 

At longer times a diffusive regime is reached, with the Einstein relation for each of the 
three Cartesian coordinates because of isotropy: 

(x 2 ) -> 2Dt; (y 2 ) -> 2Dt; (z 2 ) -> 2Dt. (8) 

We indeed find no difference between the three components of the bulk. From linear in- 



terpolation of this line in the log-log graph one reads = 0.034cr we/m, in good agreement 
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model 


source 






r* 


D* 


D* 


D* ± 


SA 


this work 


0.212 


1.17 


13.2 


0.020 


0.038 


0.04 




this work 


0.678 


0.83 


4.3 


0.034 


0.13 


0.10 


LJ 


Liu et alJ£ 


0.75+ 


0.83 




0.037 


0.15 


0.075 




Meier et alM 


0.678* 


0.83 




0.033 






LJ 


this work 


0.848 


0.74 


2.27 


0.082 


0.20 


0.16 




Meier et al4£ 


0.848* 


0.74 




0.078 







TABLE I: Table of results, in reduced units. Listed: type of model (SA: soft alkali model of 



Ref. 
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al, Ref. 



, source of the data quoted ((bulk values have been interpolated from data of Meier et 
45l ). reduced temperature T* = k^T/e, reduced density of the liquid phase p*- = pii q <r 3 , 



reduced residence time r* = Ty^e/m/a, and reduced diffusion coefficients D* = D^m/e/a: for 
the bulk (-£>*,), parallel to the liquid-vapor interface (-Dm), and perpendicular to it (D*_). 



with previous data,— >^ see Table HI 

The time for the crossover between the ballistic and the diffusive regimes can be estimated 
from the crossing of the two linear regression lines. In this case this is t c m 0.099a a/ m/e. 
The corresponding mean displacement in each direction would be Ax c ~ 0.081a if read from 
the intercept of the lines, or « 0.074a from the MSD curves at t c . We also list these numbers 
in Table HU where we will also include results for other temperatures that will be discussed 
in the next Section. Since the density is p\ iq = 0.83a~ 3 , assuming a local coordination close 
to the FCC packing, the typical intermolecular distance would be d — (v^/p) 1 ^ 3 ~ 1.2a, 
and collisions would take place with typical displacements around Ax ~ d — ao, where <Jq 
is the minimum of the LJ potential, 2 1//6 a. This prediction yields Ax ~ 0.072a, consistent 
with the values found. 



B. Diffusion parallel to the interface 

Turning to the molecules at the interface, i.e. those selected as surface molecules by the 
intrinsic sampling method, their MDS, along the three Cartesian coordinates also plotted 
in Figure HI The ballistic regime at short times is the same for all three components, and 
the same as for the bulk; but at longer times the curves are very different from the bulk, 
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the x and y, parallel, components remaining indistinguishable, the z component differing. 
We focus on the parallel diffusion in this section, for which we could expect the first two 
relations of Eq. to hold, but with a different value of the diffusion coefficient: 

(x 2 ) - 2D„t; (y 2 ) - 2D„t. (9) 

At long times, the MSD for the parallel components of our surface molecules show this 



expected linear growth with t limit, and we find D\\ = 0.13cxye/m, Therefore, the parallel 
diffusion coefficient is almost four times larger than the bulk one. This remarkable difference, 
with similar factors, had already been reported in works on water ethanol,— dimethyl 
sulfoxide,- and liquid-liquid interfaces in LJ mixtures&ii (although not all of these works 
discriminate different components of the diffusion coefficient) — we should remark that, on 
the other hand, very little change has been obtained for liquid-liquid interfaces in mixtures 
of water and other polar liquids.- 1 ^ Two works consider the liquid- vapor interface of the LJ 



fluid: a value of Dm = 0.130" ye/m, at a similar temperature of k-^T = 0.75e, is reported in 
Ref. 12 (also included in Tabled), in good agreement with our result. On the other hand, 
Ref. |20| find only a twofold increase over the bulk diffusion coefficient, even if their choice of 
temperature is again very close, k^T = 0.75e. 

The time for the crossover between the ballistic and the diffusive regimes is t c ~ 



0.38(xym/e. The corresponding mean displacement in each is either Ax c ~ 0.32<r (from 
the intercept of the lines), or « 0.23a (from the MSD curves). This means that the same 
approximate increase with respect to the bulk by a factor of about three applies to the 
characteristic time between collisions, the characteristic length between collisions, and the 
diffusion coefficient (this is of course consistent with (x 2 ) oc Dt). 

As shown in Figure HI the crossover from the ballistic to the diffusive behavior is quite 
different for the bulk and the surface molecules. The MSD in the bulk converges to the 
diffusive regime from above, i.e. after leaving the rapid ballistic regime, there is a time 
interval in which the relative growth of the MSD is slower than in the asymptotic diffu- 
sive regime. On the contrary, the parallel diffusion of the surface molecules shows a much 
smoother interpolation between the two limiting regimes, which may be interpreted as a 
signature of the stronger disorder in the correlation structure at the surface, causing a wider 
time distribution for molecular rearrangement leading to the diffusive regime. The MSD is 
always below the asymptotic diffusion, and with very little difference between the normal 
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(z), and transverse (x,y) directions up to t ~ 0.4<r a/ m/e, and displacements ~ 0.2cx, which 
are approximately half way between the ballistic prediction and the observed results for bulk 
for the same t. These values are also listed in Table HTl under the heading of "split." 

The surface molecules reach the diffusive regime only for typical displacements larger 
than one molecular diameter a, which require a time similar to the residence time. If we 
compute the time needed to diffuse to a displacement of a we find t a = 4.2crym/e, a value 
very close to that of r (4.3 in reduced units). This means that only about one third of the 
molecules at the surface (~ 0.8 x e~ x ) remain in it before diffusing a transverse distance 
similar to their diameter. We have to bear in mind the exponential decay for the number of 
molecules which remain at the surface after a time t — e.g. a molecule would likely move a 
distance 2cr on the surface in a time t ~ 4r, while in the bulk it would need typical times 
10 times larger to diffuse the same distance. However, less than two percent of the surface 
molecules (~ 0.8 x e~ 4 ) would remain as such for t = At and longer times, so that the large 
value of D\\ has limited relevance for the actual surface kinetics. The turnover process of 
molecules from the surface to the bulk phases, and its reverse, would be at least as relevant 
as the diffusion on the surface. Therefore, it is most important to analyze that turnover 
process, beyond its simple description in terms of the typical time r. In the next subsection 
we show how the intrinsic sampling method may also give an effective diffusion coefficient 
for the movement of the surface molecules in the z direction, in their wandering which will 
eventually take them out of the surface. 

C. Diffusion perpendicular to the surface 

The naive relation 

(z 2 ) -> 2D ± t (10) 

is bound to fail if we restrict the averaging to the molecules at the surface — these are 
limited in space to the region about which the IS fluctuates (both in position and in time). 
The MSD for the z coordinate will therefore tend to a constant value equal to the squared 
width of the surface molecules density profile. An intermediate diffusive regime between 
ballistic behavior and this final plateau may appear sometimes, but there is no reason to 
expect this in general. Indeed, this is not the case in the situation considered here, as is 
obvious from the dashed line in Figure HI 
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We show how to go beyond this naive prediction in a series of steps. First, Einstein's 
equation is more general than Eq. fflOl) : this would be the second moment of a probability 
distribution: 



P{z,t) 



y/4ir D±t 



cxp 



4D ± t 



(11) 



(12) 



This, of course, is the solution to the diffusion equation 

dP(z,t) _ D d*P(z,t) 
dt dz 2 

with an initial Dirac delta function distribution: P(z,t = 0) = S(z), and vanishing values 
of P and its flux for large ±z. In our MD simulations that P(z, t) would correspond to 
selecting molecules which at the time t = were within a very thin slab in the z direction, 
and which in terms of their displacement z, would diffuse to have a probability distribution 
P(z, t) after a time t. 

The same experiment could be done selecting the molecules on the intrinsic surface at 
time t = 0, and following their z displacement with t, but that would soon mix the surface 
and bulk diffusion as commented in the introduction. Alternatively, we may represent in 
P(z,t) only the z position of those molecules which have been continuously at the IS from 
t = 0. That probability distribution would not expand indefinitely, as required by the 
constant value of the MSD for the z component at long times in Figure HI That effect can 
be ascribed to an external potential U(z) that constrains surface molecules to a particular 
location in space. Indeed, this potential should be identified with a potential of mean force 
(PMF): the mean-field description of the action of all other molecules on the one diffusing 
on the surface. Our diffusion equation would now be a Smoluchowski equation: 



dP(z,t) _ D d_ 
dt dz 



dP(z,t) + ^ pM dU(z) 



(13) 



dz k^T ' dz 
(application of Smoluchowski equations to diffusion in non homogeneous media can also 



be found in Refs. 
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j28l ). The solution to this equation in the steady state would be an 
equilibrium distribution 

U(z) 



P eq (z) = P(z, t — > oo) oc exp 



(14) 



A final ingredient comes from the fact that molecules are continously leaving the IS, a 
fact which may be incorporated through a loss term in the Smoluchowski equation: 



dP(z,t) 
dt 



D — 

± dz 



dP{z,t) 
dz 



1 P M *M 



dz 



U i0 (z)P(z,t), 



(15) 
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where we have to plug the input /output frequency at each z position, v- 10 (z), discussed 
above and shown in Figure [3) The function P(z,t) solving this equation would not keep its 
normalization, since its z integral will decay with time, as the probability that an initially 
selected surface molecule would remain at the surface at time t. 

For very long times we expect that, independently of the way we initially select them, 
the remaining surface molecules will factorize as in Equation (jlj): P(z,t) ~ exp(— t/r)P (z), 
in terms of the probability distribution for the z coordinates of old surface molecules, in 
Figure [2, and decaying with the typical time r, so that (|T5|) becomes 



D — 

± dz 



dP Q (z) | P Q {z) U(z) | 
dz fccT dz 



v i0 {z)P (z) - - 

T 



0, (16) 



which, for a given D± could be used to get the effective potential of mean force which is 
compatible with forms of P {z) and fi {z) obtained from our MD simulations, and related 
to the turnover time r through (jH]). The potential of mean force could be easily obtained by 
integration of Equation (fl6l) : 



1 dU(z) d(logP ) + _^_ , <h , 



I l\ 1 

V\o{z ) 

r 



P (z>). (17) 



k^T dz dz .D_|_ 

Thus, the only unknown in Equation ( Tl5l) is the parameter D±, which is required both to 
get U(z) and to set the scale in the time variation of P(z,t). If we integrate that equation 
from the equilibrium probability distribution for all the surface molecules, it would evolve in 
time towards the exponentially decaying distribution of old surface molecules, and we could 
get D± through the comparison of the actual evolution of P(z,t) with the solution of the 
generalized Smoluchowski equation. A better estimation of D± may be obtained if we start 
with a narrower probability distribution P(z, 0), representing the surface molecules which are 
initially within a thin slab around a given zq. The evolution of this probability distribution 
shows a first stage of broadening, dominated by the diffusion so that ((z — z ) 2 ) — > 2D±t. 
In the later stage, the confining effects of U(z), and the losses regulated by u io will saturate 
the MDS and lead to the mature distribution P(z,t) ~ exp(— t/r)P (z), with a prefactor 
which would depend on the value of D±. The description of both the early and late stages 
with the same effective diffusion coefficient, which should also be fairly independent of the 
initial zq position of the chosen surface molecules, is a strong requirement to confirm the 
validity of the whole scheme, and hence the relevance of the best fitting value of as a 
true normal surface diffusion coefficient. 
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In Figure [5] we show results for the time evolution of a peak initially at zq = —a from the 
maximum of P (z) (its negative sign meaning toward the vapor side), with an initial width 
of 0.2a. The whole evolution of the P(z, t) can be traced, but we will just focus on its main 
features: its mean position, its MSD, and its decay (normalization). We try to fit these 
three curves with the numerical solution of Eq. (1151) . and we obtain a very good agreement 
fitting a perpendicular diffusion coefficient D± = (0.10 ± 0.01)o~^/e/m that is comparable 
to, but smaller than the parallel one. Of course, the exponential decay of N(t) reproduces 
the previous value of r = 4.26cr \fmje (see dotted line in Fig. [5]). We have checked that our 
results are largely independent of the initial position zq, except for peaks very close to the 
liquid phase. Indeed, the procedure seems to fail somewhat on this side, presumably due 
to the more involved nature of the process close to the liquid, with "interstitial" molecules, 
fluctuating rapidly (ballistically) from being considered as belonging or not to the IS. Those 
are again the expected limitations of any attempt to locate the outmost liquid layer, and 
we should be ready to accept some degree of "softness" in that concept. Nevertheless, the 
whole procedure provides a fairly well defined value of the normal surface diffusion, with a 
value between those of the bulk and the tangential surface diffusion, and close to the later. 

Remarkably, there is only one work, to our knowledge, that discusses normal diffusion in 
a liquid- vapor interface, Ref. [l2|. They report a value (also included in Tabled) that is lower 
than ours, and closer to the bulk value. This is presumably due to their use of slabs, by 
which some of slower bulk diffusion is mixed with the normal diffusion. The remaining works 
all discuss liquid-liquid interfaces, providing very similar values for both diffusion coefficients 
at the interface,- 1 ^ or lower values for the normal coefficient^ 1 ^ (even lower, in fact, than 
the bulk values). It is our feeling that this discrepancy stems mainly from the saturation of 
the curve corresponding to normal diffusion that has been discussed (short dashed curve in 
Figure H]): if a diffusion is computed for times shorter than the split time, a similar value for 
both coefficients will be found; but for longer times, an (effective) lower normal coefficient 
will be obtained. 



IV. DISCUSSION 

In order to support our main claims, we have performed simulations at other tempera- 
tures. We would like to consider temperatures higher and lower than the one considered. 
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FIG. 5: Diffusion of a narrow peak, comparing simulation results (solid lines) and predictions from 
Eq. (|15p (dashed lines). Upper left: mean position (first moment of the distributions), in units of 
a, lower left: MSD (second moment of the distributions), in units of <r 2 , right: fraction of surviving 
molecules N s (t)/N s (0) (normalization of the distributions). Dotted line in the later: rescaling of 
the dashed line to demonstrate we obtain the proper residence time. All graphs are versus reduced 



time t (units of (j\J m/e). 
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TABLE II: Table of relevant times and lengths for diffusion at different temperatures, in reduced 



units: ym/ecr for time, a for lengths, e/ks for temperature. Listed: model type, temperature, 
residence time (these first two columns appear also in Table H]), crossover time and displacements 
for the bulk, crossover time and displacements for the intrinsic surface, split time and displace- 
ment, time to cover a length of a. For the crossovers, two displacements are listed, one from the 
intersection of the limiting lines, another from the curves themselves. Details in the main text. 



Hence, we have collected results for the LJ liquid-vapor interface at a temperature 25% 
higher, k^T = 0.848e. It is not possible in principle to lower the temperature below the 



triple point of LJ — hence, in order to observe the influence o 
ture, we have considered the soft alkali (SA) potential of Ref. 



a sizable drop in tempera- 



35l . The parameters in this 
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model have a similar meaning to the ones in the LJ: a is the distance at which the poten- 
tial vanishes and e is related to the volume integral of the intermolecular potential function. 
Thus, the estimated critical temperature for this model is k^T c m 1.43e, similar to the one of 
the LJ fluid (k^T c « 1.33e). This is a potential engineered to result in a very low triple point 
temperature of ksT t m 0.15e. We have considered a temperature of k^,T = 0.212e. The lat- 
eral size is L = 9.025a, and the number of components employed in the IS Fourier analysis is 
8. In addition, the optimum value of the surface density is, for this model n s = 0.7a~ 2 . The 
curves showing the square root of the MSD as a function of time, for the two temperatures, 
are given in Figure EJ 

The results for the residence time and for the diffusion coefficients are included in Table [B 
We see that with increasing temperature decay diffusion becomes faster and residence times 
shorter, as is to be expected. The ratio D^/D^ is about 2.4 for the higher temperature (LJ 
fluid), since the increase in Dm is less drastic than that in D^. For the lower temperature 
(SA fluid) this ratio is about 2: in this case the reduction in D\\ turns out to be more drastic 
than that in D^. This is probably due to the fluid phase being already very dense, while 
the molecules at the surface still have the option of hopping along it — surface diffusion 
is therefore reduced due to the prevalence of energetic attraction (exerted mainly from the 
bulk liquid) over entropy at lower temperatures. 

We have also repeated the analysis of typical times and lengths for the curves in Figure 
[6l and collected this information in Table [III The most obvious correlation is between the 
time needed to diffuse a length of a, t a and the residence time. In all three cases, these two 
quantities are nearly equal. This corresponds to a natural assumption of a certain isotropy 
in the diffusion process: by the time a molecule at the IS has displaced a length of one 
molecular diameter, it is likely that the molecule has also left the IS (which is of course 
consistent with the IS being of molecular width). In other words, the product of D\\ and r 
should be nearly constant, which indeed is true, providing a typical length of y/Dp- ~ 0.7a 
in all three cases. 

There seems to be some agreement in the split length, the approximate displacement 
beyond which the parallel and perpendicular diffusion curves split, with a value around 0.2a 
in all cases. Other correlations are seen to be model-dependent. For example, the trend for 
typical bulk crossover times and lengths differs in the SA model due to its shallower potential 
well. Indeed, the argument given above in terms of a local coordination close to the FCC 
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FIG. 6: Same as Figure E] for a lower temperature of k^T = 0.212e, for the SA model (left), and 
a higher one of k^T = 0.848e for the same model, the LJ fluid (right). Solid lines: x, y, and z 
components in the bulk liquid (hardly distinguishable), long-dashed lines: x and y components 
for the intrinsic surface (hardly distinguishable, except at the longest times due to statistical 
noise), short-dashed line: z component for the intrinsic surface. Dotted lines: prediction from the 
Maxwellian velocity distribution, Eq. ([7]), and fits to Einsteinian diffusion equations, Eq. (|8|). 

packing is still valid for the higher temperature, yielding a value of Ax ~ d — <jq = 0.12<r 
(to compare against 0.14a), but fails for the SA model, with a value of d that is below the 
appropriate cr « 1.5a. 

In addition to varying the temperature, it is natural to explore the dependence of diffusion 
with the interfacial area. On one hand, the CW framework leads us to think that there 
could be a change of D± with the area. On the other hand, if this quantity is a true intrinsic 
property of the surface, it should be area-independent. In order to explore this possibility, 
we have scaled the lateral length of our simulation box, L, by factors of 1/2, 2/3, 3/2, and 
2 (i.e., the area going from 1/4 to 4 times the original). For the last two cases, we have 
increased the number of molecules to 8000 and 15625, respectively (otherwise, the liquid film 
would have been too thin); we have also enlarged the cell in the z direction, with L z = 200 
for the smallest area in order to accommodate a thicker liquid slab. In order to keep the 
same resolution in the IS Fourier treatment, the 12 functions for the original area will now 
be 6, 8, 18, and 24, respectively. 

As shown in Figure [7J functions P {z) and z/ io (z) contract for smaller areas and stretch 
for larger ones, as is to be expected from the general broadening due to CWs. (In order 
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zla z/a 



FIG. 7: Comparison of the different distribution functions for different lateral sizes. Left: in- 
put/output rate per molecule, Vi (z); right: density profile for old surface molecules, P {z). Solid 
lines: original lateral size L = 10.46<r (same curves as in Fig. dashed line: lateral size twice 
longer; dotted lines: lateral size twice shorter. All functions are normalized to unity. 

to obtain these curves, it is important that the various properties of the IS should be 
computed by subtracting the z value of the q = (constant) mode — otherwise, bulk 
density fluctuations produce a spurious broadening of the distributions at smaller areas,— 
contrary to the contraction that is found). Nevertheless, the repetition of the analysis by 
means of the Smoluchowski Equation (jT5l) . shows no measurable dependence of D± with 
the area within our 10% error bar. This supports the identification of D± as a true intrinsic 
property of the surface, and confirms the methods described. Indeed, simpler methods based 
on dividing the system in slabs,— ^ or selected "outmost" molecules^^i would have resulted 
in a clear dependence of D± with the area. 



V. CONCLUSIONS AND FUTURE WORK 



We have described an application of the intrinsic sampling method to the analysis of 
dynamical processes at the liquid-vapor interface. The main conclusion is that a liquid 
surface is a region of enhanced molecular mobility, with respect to that in the bulk liquid, 
but without strong anisotropy. The characterization of the normal diffusion constant for 
the molecules at the outmost liquid layer requires a much more elaborated method that for 
the transverse component, but at the end of the day we get similar values for D± and Dh, 
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both at high and low temperatures for the simple liquid models studies in this article. This 
fact is consistent with our results for the residence time, which governs turnover rate, i.e., 
the rate at which molecules enter and leave the intrinsic surface — a time that turns out 
to be comparable to the typical times for the threshold of diffusion. Thus, the typical time 
required for a molecule to travel a distance on the order of a molecular diameter is very 
similar to its residence time as a surface molecule. Therefore, even if diffusion coefficients 
can still be computed for molecules that stay at the surface for times long enough, the 
turnover processes are equally important when discussing the dynamical properties of the 
interface. We have discussed two main features of these processes: the overall residence 
time, and the input/output rate per molecule, Vi (z), the spacial distribution function of the 
turnover rate. 

The particular details of interfacial dynamics will be, of course, model-dependent. We 
next provide some relevant cases on which the method describe here could be applied. 
References will be given to previous works on these systems, but we would like to emphasize 
that these works have employed the usual approach (by means of slabs) , therefore the present 
approach (by means of the IS) could shed new light on the structure and dynamics of systems 
of considerable applied interest. 

The most important applications of this technique would now be realistic models of 
complex fluids, such as the liquid- vapor interface of water— ^ In this case, the more orderly 
nature of liquid phase, and the high surface tension, are likely to provide additional stability 
for the molecules at the interface, thereby leading to longer residence times and more surface 
molecules reaching the diffusion regime. We have already started an effort to study this 
system from the point of view presented in this article — in particular, the discrepancy in our 
value of the normal diffusion coefficient for the LJ fluid with results of Ref.Q suggest that the 
corresponding values for water will likewise differ. We also intend to employ this approach 
to clarify the controversy surrounding normal diffusion in liquid-liquid interfaces.- 1 *^^ 

Similar systems of interest include aqueous solutions^^^ 1 ^ (the later is a useful review 
on ion solvation at liquid surfaces) and surfactants at interfaces.— MJJLZL Surfactants will be 
usually located at the surface, and most will reach the diffusive regime. The relationship 
between structure and dynamics is specially interesting in systems with strong profiling, such 
as liquid metals, 20 confined water,— i22i2L2Ii2S and liquid-solid interfaces. 29 Similarly, diffusion 
in amphiphilic bilayer structures,—" 2 ^ micelles,— and microemulsions can be described in a 
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manner similar to the one presented here. 
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